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Abstract 

The two-equation k — ( turbulence model of 
Chien has been implemented in the WIND Navier- 
Stokes flow solver. Details of the numerical solution 
algorithm, initialization procedure, and stability 
enhancements are described. Results obtained with 
this version of the model are compared with those 
from the Chien k — € model in the NPARC Navier- 
Stokes code and from the WIND SST model for 
three validation cases: the incompressible flow over 
a smooth flat plate, the incompressible flow over 
a backward facing step, and the shock-induced 
flow separation inside a transonic diffuser. The 
k — € model results indicate that the WIND 
model functions very similarly to that in NPARC 1 , 
though the WIND code appears to be slightly more 
accurate in the treatment of the near-wall region. 
Comparisons of the k — e model results with those 
from the SST model were less definitive, as each 
model exhibited strengths and weaknesses for each 
particular case. 


Nomenclature 


a 

speed of sound 

Cf 

local skin friction coefficient 

E 

total energy per unit volume 

H 

step height, diffuser height 

ir 

diffuser throat height 

ij.k 

computational coordinate indices 

k.k+ 

turbulent kinetic energy, = Ar/u; 

Mt 

turbulent Mach number 

P.P, 

static, total pressure 

Po 

plenum pressure 

Rr 

Reynolds number based on reference 
conditions, a^.Ljv^ 
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Ret 

Reynolds number based on turbulence 
quantities 

Re r 

Reynolds number based on axial 
position 

T, T t 

static, total temperature 

t 

time 

(/, T, IT 

contravarient velocity components 

u, u, w 

Cartesian velocity components 

Mref (k — t) 

reference velocity for turbulent kinetic 

energy limiter 

U T 

friction velocity, \Jr w jp 

i/+ 

normalized velocity, u/u T 

— wu,— t/i J + Reynolds stress, —uv + — —uv/u~. 


Cartesian coordinates 

x r 

reattachment location 

y + 

normalized distance from wall. yu T ju 


rate of dissipation of turbulent 
kinetic energy, = vtju A T 

7 

ratio of specific heats 

/u/o 

laminar, turbulent viscosity 

V 

laminar kinematic viscosity, pj p 

n 

production of turbulent kinetic energy 

p 

density 


turbulent Prandtl numbers 


curvilinear coordinates 

Difference 

Operators 

ij,k 

= (Qi+\j.k-Qij.h)lU 




= (Qi+ij.*-Qwij.k)/2A* 

hQij,k 





Introduction 

In 1992, NASA Lewis Research Center (LeRC) 
and Arnold Engineering Development Center 
(AEDC) formed a partnership to enhance the 
military and commercial competitiveness of the 
Unites States through the establishment of the 
NPARC computational fluid dynamics code. This 
NPARC 1 Alliance has developed and supported the 
code by drawing from the talents of individuals at 
both centers as well as from other organizations 
outside the Alliance. 
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As explained in reference 1, recent, events have 
led the Alliance to undertake the formidable task 
of combining the predictive capabilities of three 
Navier-Stokes flow solvers into a new code called 
WIND. These solvers are the NPARC code of 
the NPARC Alliance, the NASTD code from the 
McDonnell Douglas Corporation which is now part 
of the Boeing Company, and the NXAIR code. The 
NASTD code was selected as the foundation for 
the new WIND code, primarily because it offered 
the most features of the three codes. Much of the 
Alliance's work for the past two years has focused 
on incorporating the desirable feat ures of the other 
two codes into the NASTD framework. 

At the time of the merger, the NASTD 
code offered a variety of turbulence modeling 
options including the algebraic models of P.I). 
Thomas, Cebeci-Smith, and Baldwin-Lomax, the 
one-equation models of Baldwin-Barth and Spalart- 
Allmaras, and the two-equation Shear-Stress Trans- 
port (SST) model of Menter 2 . A variety of k — t 
models had also been incorporated into the code, but 
were only moderately stable and required the user 
to be well- versed in the art of turbulence modeling. 
These k — € models were subsequently removed from 
the code prior to the merger activity. 

Meanwhile, users of the NPARC code have had 
success with the Chien 3 k — c model. The NPARC 
implementation has proven to be relatively robust 
and numerically stable for many types of flows. Of 
the low Reynolds number k — e models evaluated by 
Patel, Rodi, and Scheuerer 4 , the models of Chien, 
Launder-Sharma, and Lam-Bremhorst. were found 
to perform the best and yielded comparable results. 
For adverse pressure gradient flows, Wilcox 5 showed 
the Chien and Launder-Sharma models to be the 
best of the k — < models. For these reasons it was 
decided to include the NPARC Chien k — c model 
into the WIND code. 

The purpose of this report is threefold: (1) Com- 
bine all of tin* modifications of previous developers 
into a single complete and concise reference, 
(2) Present several validation cases for the present 
implementation. (3) Compare results from the Chien 
k — ( models in the NPARC and WIND flow solvers 
with those from the SST model in the WIND code. 

This will hopefully aid users of the NPARC" code 
in transitioning to WIND by providing a direct code 
to code comparison of the results obtained using the 
same computational mesh and turbulence model. In 
addition, results from the WIND SST model will also 
be included to demonstrate some of the strengths 
and weaknesses of each model. 


Numerical Algorithm 


Development of the Chien k-c model in the 
NPARC" flow solver has been aided by several au- 
thors over the years. The basic algorithm for solving 
the turbulent transport equations is described by 
Nichols 6 with some stability enhancements added 
by Georgiadis, Chitsomboon, and Zhu ' . 

For the WIND implementation, the k — c 
equations are nondimensionalized in a manner 
consistent with the mean-flow equations in the 
following way: 


X 

= x'/V 

p 

= P'hPL 

p 

= p’/fL 

T 

= T'/V^ 

u 

= «7a'oo 

E 

= E'/p'^a'l 

p 

= P'/p'cv 

k 

= k‘/a'l 

Pt 

= p't/p'o* 

e 



Transforming to generalized curvilinear coordi- 
nates, the k — e equations can be written as 


Q 

Fi 


(h 


S 


n 


OQ m dGj 

Or + cK, <K, + 5 


J- 1 

J- 1 


pk 

P € 

pUik 

pUit 


{J Re)' 1 


Pk'VZi ■ V* 
/'.VC ■ Vf 


(JRe)- 1 


n-pe (1 + F (M,)) Re - 2pX 
Oi/iIH - C\-,h^-Re - 2 



duj duj 
dxj dxj 


2 dui 

3 dxj 


+ b ’ jRt 


Pk = P + Pt/Vk 

p ( = p + p t /(T c 

ok 2 

Pt = C^fJ- — Re 

pc 

/, = 1.0 

h = 1 0 - 0.22 exp ( Rc t /6) 2 ^ 

/„ = 1.0 -exp (-0.011%+) 


where C t , = 0.09, C ( \ = 1.35, = 1.80. cr k . = 1.0, 

o ( = 1.3. 
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Early results with the k — c model in the NPARO 
code revealed that, nonphysical negative production 
terms could arise through the use of the complete 
production term. Thus, terms on the second line of 
the II definition are neglected, which corresponds to 
use of the incompressible form. This same approach 
has been used in the WIND code. 

To enhance predictions at higher Mach numbers, 
a compressibility correction of the form 

F (Mt) = a* • M ax (Mjf — Mf o f 0) 

has been added, and selection of the Sarkar^ or 
Wilcox y compressibility corrections is done through 
the specification of the constants a*, and M 1 0 as 
indicated in Table 1. The turbulent Mach number 
used in these corrections is defined as A if — 2 k/a 2 , 
where a is the local speed of sound. Note that 
the compressibility correction can be turned off by 
setting Qk to 0. 

Table 1: Constants used in selecting the 
com pressibi 1 i t v cor rection . 


(±k 

Mto 

Correction Type 

1.0 

0.00 

Sarkar 

1.5 

0.25 

Wilcox 

0.0 

0.00 

None 


Both the Sarkar and Wilcox compressibility 
corrections are designed to improve the prediction of 
compressible flow jets by including the compressible 
portion of the dissipation rate in the transport 
equation for the turbulent kinetic energy. These 
corrections use simple algebraic relations between 
the solenoidal and compressible dissipation rates. 
The effect of these corrections is to reduce 
the turbulent kinetic energy in compressible flow 
regions. In terms of supersonic jet predictions, this 
results in slower jet spreading rates, reduced mixing, 
and a longer core length. 

The turbulent transport equations described 
above are solved decoupled from the main flow solver 
using an approximate factorization approach 

RHS = At ( F™ - c) n F? - < 9 C F” 

+ 0 e G n x + d v GV> Hh + S" ) 


[/ + Ar (V 4 i+ + A 4 Af - ) ] A Q**= RHS 

[/ + Ar (V„i+ + A„A~ - b„B> - f’)] AQ* = A <?** 

[/ + Ar (V„A+ + A c ,47 - ^B..) ] AQ = AQ* 

where the operators operate through 

AQ**. A, B, and C are the .Jacobian matrices 
resulting from the linearization of F. (7, and S. 


Af = 
B, = 
C = 

T 7± _ { J 


\u? 

0 

0 

rf 

' Ml 

0 

0 

b\. 

’ Cll 

c 12 

. f ' 21 

C 22 

± \Bi\) /2 


bU 

b'h 




6 *„ = 


in- 


,1 Hr (£+$ + £) 6 t (j 
>H Wi + % + 9-) b>) (7 


J lit 

in- 


Jlie (C* + Cy + C;’) f- 


tL b i 

°n 

Vk 


c\i = 


ci 2 = 


C21 = 


1 

Fe 

l 


2 n _ 2 ji_ 

pk py - 


-ri 


Re [ pc 

1 


- (1 7 * F(M t ))Rc 


Re 


(C t xhW-FC €2 f 2 pcRe) — 

pk- 


-n-nA'.^Rt^)- try 


M 

36 


Rf\ 


1 


-2C a h j Re: - ( - V 

k py- \ l 


+,u2 S Cl J (~le 1 


The convection terms on the RHS are discretized 
using the Total Variation Diminishing (TV I)) 
upwinding scheme of Gorski 10 , which may be first, 
second, or third order accurate depending on the 
constants used. Because the convection terms 
contain first-order derivatives, they are similar to 
a system of linear hyperbolic equations. The use of 
a standard central-differencing scheme can result in 
nonphysical oscillations in the dependant variables, 
especially in regions of high gradients. The TVD 
property results in solutions which are essentially 
oscillation-free. 


The scheme begins by representing the convec- 
tive terms with numerical fluxes 




Ji+l/2 - /t — 1/a 


A* 
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A first-order numerical flux is then defined as 

*«'+l/2 = 2 (/* + ! + /' - d fi + 1/2 + 

Second- and third-order flaxes are obtained by 
adding corrections to the first-order flux. 

f‘+ 1/2 = /, ' + l/2 - d ff+3/2 ~ d fi+ 1/2 

+ ^<+1/2 + V^-,/2 

The symbols df and d/ denote flux limited values 
of df and are computed using the minmod operators 
described by Chakravarthy 11 

d ff+ 3/2 = minmod (df~ +3/ .,, 3 ■ df ~ +1/2 ) 

f/ /r+i/2 = (d-i/2’-^ ‘ Ch3/->) 

^/h-1/2 = minmod \dff +l/2 ,.3 ■ dff_ , /2 ) 

d ft 1/2 ~ minmod (rf/^. 1/2 , 3 • df ? + 1/2 ) 

■mm mod ( <t\ ,t/) = shyn (j*) • 

T7j.«jr {0. m/n [|j*| , y ■ ( ar )] } 

where the compression parameter, /i, is defined as 
3 — 0 


where 


0i+l/2,j,Jfe 

= [ 




= [ 

O' (Cr kx “1" C'/ 

ky T 

C: k t )] ; + , /2 



1 ( 

pt v 


r *t + l/2,j,fc 


_J Re V' + 

* -)\ 

1+1/2, J,* 

( )j + l f 2 .j,k 

= (&+{ + rixk n 

T Ch 

H‘ h+l/ 2 , j,k 

(ky) i + l f 2 j k 

— (£, 1 J k/ “f- Ijykff 

+ C,v k'C )i + ]/ 2 ,j'k 

)* + l/2,j ,k 

= (Zzk ( + r) t k n 

T C ^ 

‘ C )? + 

and // and p t 

at i- j- 1/2 are computed from a si 


average between / and i T 1. To determine the 
gradient of k in the sweep direction, 

and for the non-sweep directions, 

2 [ ( M + (M,+i] , fc 


( ^ V ) i + 1 / 2 ,j , fc 
(*C )| + l/2,i 1 Jb 


IK +(*<>*,]., 




/J = 


1 “ 0 


~ 2 (^ A ’ + 1 1 )*,j 


Second- or third-order schemes can be obtained 
by setting o equal to -1. or 1/3 respectively. Then, 
for example, the convective terms in the turbulent 
kinetic energy equation can be written as, 

/ = a k 

a = J~ { pU 


The convection terms on the LHS are discretized 
using a flux-splitting technique and the diffusion 
terms are discretized using a second order central 
difference. By neglecting the cross-derivative terms 
that would normally appear in Bj , the resulting 
equations form a block tridiagonal system which can 
be readily inverted. 


and the flux differences are 


dff + j /2 = max (tt f+ |,0) • ki+\ - max (a,-, 0) • 
<lf~ +l j 2 = min (a, + i , 0) fc i+1 - min (Oj , 0) • k t 

For the convective terms in the tj— or in- 
directions, l is replaced by the appropriate 
contravarient velocity and i is replaced by the 
corresponding coordinate index. 

Since the diffusion terms are composed of second- 
order derivatives, which tend to have a smoothing 
effect, oscillations in the dependent flow variables 
is not of concern and a standard central difference 
discretization can be used. For example, the 
diffusion of turbulent kinetic energy in the £ — 
direction is computed from 


<ku 


(fh+i /2 “ <U-\r>)j k 

A* 


Initialization 

The turbulent transport variables for the Chien 
k — e model (fc,e,p*) can be initialized using two 
different techniques. The preferred method uses 
an assumption of turbulent equilibrium, in which 
the production of turbulent kinetic energy equals 
the rate of dissipation, together with an existing 
turbulent viscosity profile to initialize the k and ( 
variables. 


pk 


U/Rr 


/ 


PWt 

( (if n R f 


In order to use this technique, the code must 
first be run a few thousand iterations using another 
eddy- viscosity turbulence model. Initializing from 
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an existing turbulent viscosity profile rather than 
uniform values aids somewhat in convergence and 
improves the stability of the model by reducing 
the severe changes in turbulence values that occur 
during the first, few iterations after initialization. 

The second method initializes the turbulence 
variables to uniform values (k,c. fit) within each zone 
using the local density. This technique has not been 
found to be very robust. 

Stability Enhancements 
Relaxation 

Updated values of k< e, and fi t are relaxed for 
a set number of iterations following initialization. 
Relaxation of these variables reduces the amount 
they may change during any single iteration. 
Immediately after initialization, the allowed changes 
are significantly reduced. This restriction is then 
gradually lifted as the last relaxation iteration is 
approached. 


Limiters 

The k — c model also uses limiters within the 
interior of each zone to increase convergence and 
stability by capping the values of the turbulence 
quantities at both the high and low extremes. This is 
usually only necessary during the first few iterations 
after initialization, when the fluctuations in k and c 
tend to be the most severe. 

Nondimensional values of the minimum limiters 
have been preset to relatively small numbers. The 
values of the maximum limiters are determined by 
user inputs for the maximum allowable turbulent 
viscosity fit ma x and a turbulent reference velocity, 
from which the maximum allowable value for the 
turbulent kinetic energy is computed using: 


The maximum dissipation rate is computed from the 
turbulent viscosity relation. 

The use of these limiters can be summarized 
as follows: (1) If either k or ( falls below preset 
minimum values then both are reset to these values. 
This typically occurs in the freest ream. (2) If the 
turbulent kinetic energy exceeds k mar , then it is 
capped at this value. The dissipation rate is taken 
to be the larger of the current dissipation rate or 
c max • (3) If the turbulent viscosity exceeds fit max, 
then it is capped at this value and the turbulent 
kinetic energy is recomputed from the turbulent 


viscosity relation. The turbulent dissipation rate is 
left unchanged. 

These maximum limiters are meant, to keep the 
solution from diverging during the initial iterations 
of the solution and care must be taken that these 
limiters are not constraining the solution upon final 
convergence. Verifying that the maximum value of 
the turbulent viscosity in the flowfield is less than 
that specified for fit mar is usually sufficient,. At 
the conclusion of a set of iterations, the WIND 
code provides the user with a warning message if 
the solution is being constrained by the maximum 
limiters. 

Variable (\ 

It is well known that the baseline k — e model 
is poorly suited to adverse pressure gradient flows. 
Rodi and Scheuerer 12 demonstrated that for these 
types of flows, the rate of dissipation near solid 
boundaries is too small relative to the rate of 
production of t urbulent kinetic energy. This causes 
the model to overpredict skin friction and predict, 
flows to be attached when experimental results show 
them to be separated. The variable formulation, 
which is derived from algebraic stress modeling, is 
designed to help remedy this problem by reducing 
the turbulent viscosity in regions of the flowfield 
where the production of turbulent kinetic energy is 
significantly larger than the rate of dissipation. The 
specific formulation used is taken from Rodi l3 : 


Cp = min 0.09 


0. 10738 (0.64286 + 0. 196077?) 
[1 + 0.357 (fl-1)] 2 


As the ratio R of production to dissipation increases 
above one, the coefficient is reduced from its 
normal value of 0.09 to limit the turbulent, viscosity. 

The variable option also provides added 
stability to the k - e model , such as in the case of 
an airfoil, where the sudden deceleration of the flow 
near the leading edge would otherwise result in an 
unrealistically high rate of production. In regions of 
the flow where the turbulence is in equilibrium, i.e. 
where the production and dissipation are balanced, 
the turbulent viscosity remains unchanged. 


Numerical Results 

The focus of this validation effort, will be on 
wall-bounded flows. For a similar code to code 
comparison of mixing results for supersonic exhaust 
nozzles, the reader is referred to reference 14. 
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Flat Plate 

The incompressible flow over a smooth flat plate 
was used as an initial validation case. The flow 
being modeled is that reported by Wieghardt 15 and 
later included in the 1968 AFOSR-IFP Stanford 
Conference u \ The freest, ream Mach number in the 
simulations was set to 0.20, slightly higher than 
that in the experiment, in order to accelerate the 
convergence rate. 

Figure 1 depicts the computational domain used 
to model this flow. A ('artesian mesh with 111 
points in the axial direction and 81 points normal 
to the viscous wall was used. The first 14 grid 
points upstream of the leading edge of the plate 
were treated as an inviscid wall to provide a uniform 
profile at the leading edge location. The grid was 
packed in the streamwise direction to resolve flow 
gradients near the leading edge of the plate and 
normal to the surface to resolve the boundary layer. 
Calculations were made on a series of grids having 
t/ + values of 1, 2. 5. 10, and 30 at the first point 
off the wall. Figure 2 indicates that these values are 
representative of the maximum along the plate 
as computed from the WIND k — e solutions. These 
grids are the same ones used in the Chien k — c grid 
sensitivity study of reference 17 with the NPARC 
code. 

Figure 3 shows the computed skin friction using 
both codes. The grid sensitivity studies shown in 
Figures 3a and 3b for the NPARC and WIND k — e 
models indicate that grid independence for both 
codes is obtained using the = 2 grid. Figure 3c 
shows this to be the case for the WIND SST model 
as well. While the WIND skin friction results for 
larger values appear to be much less sensitive 
than the NPARC results, analysis of the turbulence 
quantities in the near-wall region reveal that the 
model predictions begin to break down severely as 
the mesh spacing exceeds y+ = 5. Figure 3d is a 
direct comparison of the skin friction results between 
codes using the t/ + = 1 grid. As can he seen, 
the WIND k - f results are improved and compare 
well with the WIND results using the SST model. 
The comparison of the computed velocity profiles at 
several axial stations given in Figure 4 also shows 
that the WIND code matches the experimental data 
quite well. 

Fxamiuatiou of the turbulence quantities at the 
last velocity station (Re r — 1.03*10') is given 
in Figures 5-7. Figure 5 compares the turbulent 
kinetic energy with the “average” experimental 
data assembled by Patel, Rodi, and Scheuerer 4 . 
While neither code matches this data exactly, the 


WIND results more closely match those obtained by 
Patel using the Chien model in a two-dimensional 
boundary layer code. The WTND SST model fails to 
capture the peak in the turbulent kinetic energy due 
to the form of the k — uo model used in the near- wall 
region. From the turbulent dissipation rate shown 
in Figure 6, one can see that the WIND results 
approach the correct limiting value away from the 
wall. Comparison of the turbulent shear stress in 
Figure 7 shows that only the WIND results approach 
the correct asymptotic value away from the wall. 

Backstep 

The second validation case is the incompressible 
flow over a backward facing step. As shown in 
Figure 8, this geometry has a step height to tunnel 
exit height ratio of 1:9 which helps to minimize 
the freestream pressure gradient due to sudden 
expansion. The experimental configuration of Driver 
and Seeginiller 18 also had a step height to tunnel 
width ratio of 1:12 to minimize three-dimensional 
effects. A variety of experimental measurements are 
available, including skin friction, pressure, turbulent 
normal and shear stresses, and velocity profiles 
downstream of the step. 

A 238 x 185 single-zone mesh was generated to 
model the region from x/H =-105 to +50. The grid 
was packed to the solid surfaces such that t/ + ss 1. 
Downstream of the step, 55 points were used in 
the recirculation region with ten of them placed 
within y+ & 30. The grid was also clustered in the 
streamwise direction near the recirculation region to 
improve resolution. 

Figure 9 shows the velocity profiles at several 
axial locations. Upstream of the step, all of 
the solutions are virtually identical. Within the 
recirculation region, however, there are noticable 
differences. The NPARC and WTND k - c solutions 
are nearly indistinguishible and appear to provide 
the best match to the experimental data. Use of 
the variable C\ L option in the WTND k — c mode] 
causes the flow to reattach further downstream, but 
does not predict the rest of the velocity profile as 
well as the standard k — e model. The WTND SS I 
model shows the greatest disparity compared to the 
data and predicts the flow to reattach even further 
downstream than the k — c model with the variable 
Op option. 

These findings relative to the extent of the 
recirculation region are reiterated in Figure 10, 
which shows the predicted skin friction coefficient. 
Both the NPARC and WTND k — ( models predict 
the reattachment to occur too far upstream and 


0 
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display under- and overshoots relative to the data. 
According to Avva, Smith and Singhal 19 this 
overshoot can be reduced by increasing the number 
of points below % 30. For a mesh with ten 
points inside this region, the present results agree 
with those presented by Avva. The variable 
option tends to reduce the turbulent viscosity within 
the separation region, thus making the flow appear 
more laminar-like and reducing the magnitude of the 
skin friction. The predicted reattachment location 
is also shown to move downstream. Unlike for 
the velocity profdes, the WIND SST model seems 
to provide relatively good agreement with the skin 
friction data. Table 2 lists the reattachment location 
predicted by each model. 

Table 2: Predicted Reattachment Locations 
for the Backward-Facing Step. 


Model 

Tr/H 

NPARC k - e 

5.31 

WIND k - e 

5.30 

WIND k - t Var. C » 

5.55 

WIND SST 

6.43 

Driver Experiment 

6.26 


Figure 11 shows the turbulent kinetic energy 
profiles at several axial locations. Upstream of 
the backstep, the flow is similar to that of the 
flat plate and one can again see that the peak 
value in turbulent kinetic energy is underpredicted 
by the SST model due to the form of the k - 
model used in the near-wall region. This difference 
appears to propagate downstream as the SST model 
consisently predicts a lower peak value than the 
k — e model at each axial station. As with the 
velocity profiles of Figure 9, there is close agreement 
between the NPAR.C and WIND k — e solutions. One 
can also notice the reduction in turbulent kinetic 
energy, especially within the recirculation region, 
caused by the use of the variable C ^ option. This 
occurs because the variable option effectively 
increases the turbulent dissipation rate within the 
recirculation zone. 

The corresponding Reynolds stress profiles are 
shown in Figure 12. Here again there is excellent 
agreement between the NPARC and WIND k — e 
solutions, except near x/H=0 where the WIND 
solution predicts a more rapid increase in the 
downward component of velocity. Far downstream 
the k — e model overpredicts the Reynolds stress 
which corresponds with the overprediction in skin 
friction shown in Figure 10. One can also observe 
how the variable (\ t correction noticeably reduces 
the peak turbulent viscosity (and consequently the 
Reynolds stress) within the recirculation region. 


The ability of the SST model to match the 
Reynolds stress data so well downstream of the 
reattachment location is interesting, considering the 
underprediction of velocity and turbulent kinetic 
energy in this same region. 

Transonic Diffuser 

The Sajben 20, 21 diffuser strong-shock case was 
selected as the next validation case. Figure 13 
is a schematic of the two-dimensional diffuser 
geometry. This configuration had an entrance 
to throat area ratio of 1.4, an exit to throat 
area ratio of 1.5, and a sidewall spacing of 
approximately four throat heights. Suction slots 
were placed on the side walls of the constant area 
sections upstream and downstream of the diffuser 
to minimize three-dimensional effects. Streamwise 
slots were also placed along the top corners of 
the diffuser to maintain a two-dimensional flow. 
Time-averaged static pressure distributions on the 
top and bottom diffuser walls were measured 
using pressure transducers, and separation and 
reattachment locations were obtained through the 
use of oil-flow techniques. Velocity profiles were 
obtained using a laser Doppler velocimeter. 

Although this geometry was tested both with 
and without externally applied oscillations, only the 
steady-state flow of the unexcited cases was modeled 
numerically. These flows were characterized by the 
ratio of exit static to inflow total pressure. For the 
strong-shock case this ratio was 0.72. 

This case was computed using ail 81 x 51 grid, 
which corresponds to the coarse mesh used in the 
investigation conducted by Georgiadis, Drummond 
and Leonard 22 with the PARC. 1 code. They found 
the grid to be sufficiently clustered in the vertical 
direction such that the first point off the wall resided 
inside the laminar sublayer. 

Figure 14 shows that without any turbulence 
model corrections, the WIND k — e model predicts 
the shock too far downstream and poorly matches 
the pressure distributions downstream of the 
shock. Use of the Sarkar compressibility correction 
improves the prediction of the shock location 
somewhat, but does not improve the flowfield 
solution in the downstream region. The variable 
Cp option has the same effect on the shock location 
as the Sarkar correction, but also improves the 
downstream pressure predictions. By using the 
Sarkar correction in conjunction with the variable 
C' fl option, the shock location is reasonably well 
matched and the downstream pressure distribution 
is improved. 
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Figure 15 illustrates the effect of the k — c 
model correction factors on the predicted velocity 
profiles downstream of the shock. Without any 
corrections, the model is unable to match the 
experimental data near the lower wall and the core 
flow velocity is underpredicted. Use of either the 
Sarkar compressibility correction or the variable 
option yields improved results both near the walls 
and in the core. Using both correction factors yields 
some additional improvement near the lower wall, 
while the upper wall region appears to be over- 
corrected. 

Cross-code comparisons of the strong shock 
results are made in Figures 16 and 17. Two sets 
of WIND k — e results are plotted. The first set uses 
the same correction factors (Sarkar compressibility 
correction used, but not the variable C ^ option) as 
the NPARC code and should be used to compare 
the new and old k — e implementations. The second 
series was computed using both correction factors 
to demonstrate the benefit of using the variable 
Cft option. Figure 16 shows that the WIND and 
NPARC 1 k — t model pressure distributions are 
again in close agreement, with the WIND code 
still predicting the shock location one or two grid 
points further downstream. With the addition of 
the variable C' /4 option, the shock location predicted 
using WIND agrees well with that using NPARC 
and the downstream pressure distribution compares 
better with the experimental data. The WIND SST 
model also provides improved downstream pressure 
distributions, but predicts the shock to occur further 
upstream than indicated by the experimental data. 

Conclusion 

The two-equation Chien k — € turbulence model 
has been successfully implemented in the WIND 
Navier-Stokes flow solver. Details of the numerical 
algorithm have been presented including the ini- 
tialization procedure, stability enhancements, com- 
pressibility corrections, and variable formulation. 
Results for the wall-bounded flows investigated 
herein indicate that the current implementation 
functions very similarly to that in the NPARC code, 
though the W IND code appears to be slightly more 
accurate in predicting skin friction at the wall. 
( omparison of the Chien model results with those 
from the SST model for these cases yielded no 
obvious favorite. 
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Figure 1 : Schematic of the Flat Plate Test Case. 
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Figure 2: Near-Wall Grid Spacing for Flat Plate Case. 




Figure 3a: NPARC Chien k-e Grid Sensitivity Study. Figure 3b: WIND Chien k-e Grid Sensitivity Study. 




Figure 3c: WIND SST Grid Sensitivity Study. Figure 3d: Comparison of Grid Independent Solutions. 

Figure 3: Local Skin Friction Along Flat Plate. 
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Figure 4: Velocity Profiles Along Flat Plate Using y =1 Grid 
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Figure 10: Local Skin Friction Coefficient Downstream of Backstep. 




Figure 1 1: Turbulent Kinetic Energy Profiles for Backward-Facing Step. 
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Figure 14: Effect of WIND k-e Model Corrections on Predicted Pressure Distributions 

for the Transonic Diffuser Case. 
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Figure 15: Effect of WIND k-e Model Corrections on Predicted Velocity Profiles 
for the Transonic Diffuser Case. 
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Figure 16: Pressure Distributions Along Top and Bottom Walls 
for the Transonic Diffuser Case. 
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Figure 17: Velocity Profiles at Four Axial Locations for the Transonic Diffuser Case. 


NAS A/TM— 1 999-209080 


19 



REPORT DOCUMENTATION PAGE 


Form Approved 
OMBNo. 0704-0188 


Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources 
gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 
collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson 
Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget. Paperwork Reduction Project {0704-0188), Washington, DC 20503. 


1. AGENCY USE ONLY ( Leave blank ) 


2. REPORT DATE 


April 1999 


3. REPORT TYPE AND DATES COVERED 

Technical Memorandum 


4. TITLE AND SUBTITLE 


Implementation and Validation of the Chien k-e Turbulence Model in the Wind 
Navier-Stokes Code 


AUTHOR(S) 

Dennis A. Yoder and Nicholas J. Georgiadis 


5. FUNDING NUMBERS 


WU -537-05-2 1 -00 


7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
John H. Glenn Research Center at Lewis Field 
Cleveland, Ohio 44135-3191 


8. PERFORMING ORGANIZATION 
REPORT NUMBER 


E- 11662 


9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
Washington, DC 20546-0001 


10. SPONSORING/MONITORING 
AGENCY REPORT NUMBER 

NASA TM — 1999-209080 
AIAA-99-0745 


11. SUPPLEMENTARY NOTES 

Prepared for the 37th Aerospace Sciences Meeting & Exhibit sponsored by the American Institute of Aeronautics and 
Astronautics, Reno, Nevada, January 1 1-14. 1999. Responsible person, Dennis A. Yoder, organization code 5860, (216) 
433-8716. 


12a. DISTRIBUTION/AVAILABILITY STATEMENT 

Unclassified - Unlimited 
Subject Categories: 34, 02, and 61 


12b. DISTRIBUTION CODE 


Distribution: Nonstandard 


This publication is available from the NASA Center for AeroSpace Information, (301 ) 621*0390. 


13. ABSTRACT (Maximum 200 words) 

The two equation k-e turbulence model of Chien has been implemented in the WIND Navier-Stokes flow solver. Details 
of the numerical solution algorithm, initialization procedure, and stability enhancements are described. Results obtained 
with this version of the model are compared with those from the Chien k-e model in the NPARC Navier-Stokes code and 
trom the WIND SST model for three validation cases: the incompressible flow over a smooth flat plate, the incompress- 
ible flow over a backward facing step, and the shock-induced flow separation inside a transonic diffuser. The k-e model 
results indicate that the WIND model functions very similarly to that in NPARC, though the WIND code appears to be 
slightly more accurate in the treatment of the near- wall region. Comparisons of the k-e model results with those from the 
SST model were less definitive, as each model exhibited strengths and weaknesses for each particular case. 


14. SUBJECT TERMS 


Turbulence models; Turbulent boundary layer 


15. NUMBER OF PAGES 

25 


16. PRICE CODE 


17. SECURITY CLASSIFICATION 
OF REPORT 

Unclassified 


18. SECURITY CLASSIFICATION 
OF THIS PAGE 

Unclassified 


19. SECURITY CLASSIFICATION 
OF ABSTRACT 

Unclassified 


20. LIMITATION OF ABSTRACT 


Standard Form 298 (Rev. 2-89) 

Prescribed by ANSI Std. Z39-18 
298-102 















